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ABSTRACT 

Previous theoretical work has speculated about the existence of double- 
diffusive magnetic buoyancy instabilities of a dynamically-evolving horizontal 
magnetic layer generated by the interaction of forced vertically-sheared velocity 
and a background vertical magnetic field. Here, we confirm numerically that 
if the ratio of the magnetic to thermal diffusivities is sufficiently low then such 
instabilities can indeed exist, even for high Richardson number shear flows. Mag- 
netic buoyancy may therefore occur via this mechanism for parameters that are 
likely to be relevant to the solar tachocline, where regular magnetic buoyancy 
instabilities are unlikely. 

Subject headings: hydrodynamics — MHD — Sun: magnetic fields 
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Introduction 



The generation of magnetic field by velocity shear and the field's subsequent evolution 
are of great importance to an understanding of the operation of the solar dynamo. W hile the 
curren t dynamo parad igm contains many complex interacting components (see e.g. iParker 
( 119931 ); iRempell (120061 )). an integral part of all current large-scale solar dynamo models is 
the creation of strong toroidal magnetic structures in the tachocline, the Sun's thin region 
of strong radial differential rotation separating the latitudinally differentially rotating con- 
vection zone and the solid-body-rotating radiative interior. Strong toroidal magnetic field 
is thought to be induced by the stretching action of the different ial rotation on a n y bac k- 
ground poloidal field (the fi-effect of mean-fiel d dynamo th eory; ISteenbeck et al.l (119661 )). 
Subsequently, magnetic buoyancy instabilities (jParkerl Il955l ) of the generated field are in- 
voked as the mechanism for the creation of distinct magnetic structures and their subsequent 
rise toward eventual emergence at the solar surface as active regions. 



In recent papers, IVasil fc Brummelll (120081 . 120091 ) (hereinafter VB1, VB2) show that, 
when naively using the commonly accepted paradigm for the f2-effect, it is surprisingly diffi- 
cult to initiate magnetic buoyancy instabilities if thermodynamic adjustments are assumed 
to be adiabatic. Specifically, they concluded from analytic (VB2) and numeric (VB1) calcu- 
lations that for magnetic buoyancy instabilities to operate in a thin tachocline, the velocity 
shear flow imposed to drive the system must be necessarily hydrodynamically unstable. 
Roughly speaking, VB2 found that Ri ~ Az/H p is required for magnetic buoyancy, where 
Ri is the Richardson number, Az is the vertical shear width and H p is the local pressure scale 
height. For the solar tachocline, the estimated Richardson n umber is very large (Ri ~ 10 3 - 
10 5 ) but the region is thin so that Az/H p < 1 (IGoughl 120071 ) and so the above condition for 
instability is therefore unlikely to be satisfied. The simulations of VB1 confirmed numeri- 
cally that hydrodynamically unstable imposed shear flows could generate magnetic buoyancy 
instabilities in a thin shear region, but stable shear flows did not. Energetically, this can be 
thought of as the following conundrum; for magnetic buoyancy to occur, the shear flow must 
transfer enough energy into a toroidal magnetic field for it to overcome the stable background 
stratification. If the shear can only build (through an Alfvenic process) magnetic field to 
the level of equipartition with the flow, and the shear is constrained by the stratification 
(since it is hydrodynamically stable), then it is difficult for the shear- induced magnetic field 
to overcome the constraints of the stratification. 

It is the case, however, that near the bottom of the convection zone the ratio between 
the magnetic and thermal dif fusivities, ( = is very small. It has long been recognized 
(lGilmanlll970l ; lAchesonlll979l ) that, in such circumstances, instabilities can occur that rely 
on the much greater diffusion rate of the stabilizing thermal component (see iHughed (120071 ) 
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for further discussions). Further, VB2 noted that instability might be enhanced by double- 
diffusive effects. Assuming isothermal (as opposed to adiabatic) adjustments VB2 obtained 
the much less stringent requirement that (Ri ~ Az/H p . It is therefore possible that a 
more-solar-like (high Ri) shear-generated magnetic field can become buoyantly unstable in 
a double-diffusive manner when (<cl. 

In this letter, we consider this possibility in a convectively stable atmosphere. We show 
that if ( is sufficiently small, and the background Alfvenic timescales are sufficiently slow, 
then double-diffusive instabilities can exist for parameters that are relevant to the solar 
interior. 



2. Equations and parameters 



We consider a Cartesian domain (x, y, z), where z is depth, x is the toroidal/zonal direc- 
tion, and (y, z) are the poloidal directions. Our system consists of an initially vertical uniform 
magnetic field, B z $, permeating a stratified layer of compressible fluid under the influence of 
a forcing designed to generate a target shea r flow u = (Uq(z), 0, 0). This system is similar to 
that described in VB1 and in ISilvers et al\ ( 120091 ). We solve standard non-dimensionalized 
equations for forced compressible magnetohydrodynamics governing the evolution of velocity 
u = (u,v,w), magnetic field B = (B x , B y , B z ), density p, temperature T, and pressure P 
based around a polytropic atmosphere, T (z) = 1 + 9z, po(z) = T (z) m , Pq{z) = T (z) m+1 . 
The non-dimensionalization uses T*, the temperature at the top of the domain, p* which 
is proportional to the total mass within the domain, = (c p — c v )T*p*, the fiducial pres- 
sure (given the specific heat capacities, c p ,c v ), B z0 , the imposed background magnetic field 
strength, the layer depth d, and time units of = d p l J 2 / P^ 2 . We impose stress-free velocity 
and vertical magnetic field boundary conditions at z — 0, 1 together with T(z = 0) = I and 
d z T(z = 1) = 9. The horizontal directions are periodic. 

The sys tem is solved i n a d omain with aspect ratio 2:1:1 in a similar manner to that in 
VBI and in Silvers et al. J2009h at a resolution of (roughly) 256 x 256 x 512. This aspect 
ratio is specified to allow for certain expected dynamics, e.g. we anticipate three-dimensional 
modes to have smaller scales in the y and z directions than in the x direction. The thinness 
of the forced shear region compared to the local scale height is the geometrical factor that 
identifies the model with the tachocline. The kinetic and magnetic Reynolds numbers for 
our calculations are roughly Re ~ 2000, Reu ~ 1000, respectively (ba sed on the forced shear 
magnitude and width). Estimates of the numerical degrees of freedom (jDavidsonl|2004j ) imply 
that our resolution is reasonable. However, we further satisfied ourselves that the results 
were consistent at varying resolutions. 
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The important parameters in the problem are the dimensionless thermal diffusivity, 
Ck = Kt*/ p*c p d 2 , the Prandtl number a = pc p /K (/i representing dynamic viscosity), 
the inverse Roberts number ( = rjc p p*/K — r)/n, and a = B\ T 2 /p p*d 2 which gives a 
measure of the Alfven speed along the background field in terms of the fundamental acoustic 
yelocity scale. Note th at ( is the critical parameter governing double- diffusive instability 
(faughes fc Weisslll995h . 



We add a forcing term, F = —aCKd 2 U x, in the x-momentum equation that would 
(in the absence of magnetic effects and instabilities) maintain a desired target velocity u = 
(Uo(z),0,0) against viscous decay We choose 
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to mimic the smooth radial shear transition believed to exist in the solar tachocline. The 
width is chosen sufficiently narrow that OJJq ~ (to within numerical precision) at the 
boundaries. 



An important derived parameter is the Richardson number, 
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where N (z) = \—g d z In {Pq^z) 1 ^ / po(z))} 1 / 2 is the local Brunt- Vaisala frequency of the stable 
background atmosphere, d z lIo(z) is the local turnover rate of the background shear and 
7 = c p /c v = 5/3. Ri measures the relative tendency of a shear flow to overturn fluid 
vertically compared to gravity's tendency to restore it to its original position. A large 
Ri implies that gravity is strongly stabilizin g whereas a small value can mean that shear 
instabilities are possible (IDrazin fc Reidll2004l ). The goal in this work is to obtain a magnetic 
buoyancy instability at high Richardson number, a result not found in VB1, but here we 
investigate the low ( regime where the effects of thermal stability are severely reduced. 

While the Richardson number provides a useful rough measure of the stability properties 
of our system, we also define the respective thermal and magnetic Rayleigh numbers 



Rt(z) 
Rb{z) 



No(z) 



is(z)k(z) 
gaB x (zf 
z/(z)k(z)Po(z) dz 



d 



In B x (z), 



(3) 
(4) 



w here n{z) = Ck/po(z), a nd v{z) = an(z). These definitions are consistent with those given 



in 



Hughes fc Weissl (119951 ) except that here the magnetic Rayleigh number, Rb, is given in a 
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form that anticipates three-dimensional instabilities (jNewcombl 1 1 96 ll ) . For the direct (non- 
oscillatory) type of instability, of interest here, both Rayleigh numbers are negative, and 
this corresponds to a thermal stratification that is stabilizing and magnetic gradient that is 
destabilizing. A Rayleigh number based on the total density stratification would correspond 
to .Rtotai = Rt — Rb- The key result of VB2, derived under the adiabatic thermal adjustment 
assumption, is equivalent to saying that -Rtotai is always negative (ostensibly stable) unless 
Ri <C 1. In the double- diffusive context, where thermal adjustments are closer to isothermal, 
the critical parameter is actually i?DD = Rt~C 1 Rbi where £ <C 1. It is -Rdd that is relevant 
to the stability of our system and it can potentially become positive even if Ri ^> 1. This is 
equivalent to the isothermal result derived in VB2 that required (Ri ~ Az/H, 



For the basic thermal parameters (see Table [T]), we take the polytropic index m = 
1.6 that enforces a convectively stable background polytropic stratification, and a non- 
dimensional lower-boundary heat flux of 9 = 5. With the chosen shear flow, this produced a 
minimum Richardson number of Ri = 2.96. While this value is significantly lower than we 
expect for the solar tachocline, it is high enough to guarantee that Uq(z) hydrodynamically 
stable. This fact was checked numerically with a purely hydrodynamic run. 

In order to see if double-diffusive instabilities can exist, we ran two simulations that differ 
only in the thermal diffusivity k = K/c p p* (varied through its non-dimensional counterpart 
Ck)- Thus ( and a are varied commensurately to maintain fixed viscous and magnetic 
diffusivities fi = aC K and r) = (C K . We keep aC K = 2.5 x 10~ 6 and (C K = 5.0 x 10~ 6 
for all our cases. For our first case, CI, anticipating instability, we choose ( = 5.0 x 10 -4 , 
a = 2.5 x 10~ 4 and Ck = 0.01. In our second case, C2, anticipating stability, we choose 
( = 0.01, a = 5.0 x 10~ 3 and C K = 5.0 x 10~ 4 . Thus, C in CI is 20-fold smaller than in C2. 

The parameter a = <j(C k Q (where Q is the Chandrasekhar number measuring the ratio 
of magnetic and diffusive timescales) sets the background vertical magnetic field strength. 
For CI and C2 we set a = 1.25 x 10~ 5 . We desire a 1 so that Alfven timescales are slow 
compared to the acoustic timescale, and Q ^> 1 so that Alfven timescales are fast compared 
to diffusion. CI and C2 have Q — 1.0 x 10 6 . It turns out that the Alfven timescale of 
the background magnetic field has important consequences for the dynamic stability of our 
evolving MHD configuration. To investigate this aspect, we therefore present a third case, 
C3, that is identical to CI except that a = 5.0 x 10 -5 , (Q = 4.0 x 10 6 ). C3 therefore has 
two-fold faster Alfven timescales. 
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3. Results 

In all three simulations, we begin with a purely vertical magnetic field, and allow the 
induction of a toroidal field layer by the action of the imposed velocity forcing. In the absence 
of two- and three-dimensional effects, this induction process is governed by a one-dimensional 
(mean) set of MHD equations: 

po(z)d t u = ad z B x + aC K d 2 z u - aC K d 2 z U (5) 
d t B x = d z u + (C K d 2 z B x . (6) 

Strictly speaking, the density profile in Equation (jSJ) will e yolve along w i th th e shear flow 



and magnetic field. However, unlike the work in VB1-2 and lSilvers et al\ (120091 ). the target 
shear flow here, Uq(z), has a large Richardson number. To an excellent approximation, the 
background density maintains its initial polytropic profile p(t, z) = po(z) in this initial phase. 
The evolution, as shown in Figure HJ of the mean toroidal magnetic field is therefore essen- 
tially independent of the thermal properties of the system, unless two- or three-dimensional 
instabilities occur. 

As time progresses, the peak field grows in the region of maximum forced velocity shear, 
and the field gradients strengthen. The ultimate question is whether a weak shear (Ri > 1) 
can induce strong enough gradients for a magnetic buoyancy instability to occur. Below, 
we show that two initially identically evolving one-dimensional MHD configurations can 
have completely different stability properties, based solely on the magnitude of the thermal 
diffusivity. 

In CI, the system does indeed initially evolve according to Equations (|SJ) and 
However, since ( = 5.0 x 10 -4 is small enough, the induction of mean toroidal field only 
proceeds for a finite time before an instability becomes apparent at t ~ 88. Figure (T5]) shows 
volume-rendered images at this time of the vertical velocity w(x,y,z) and the fluctuating 
toroidal magnetic field B x (x,y,z) — B x (z). In the initial induction phase (governed by 
Equations (j3J) and ©) both these quantities are zero. However, Figure (jSj) clearly shows 
that a wave-like perturbation in the vertical velocity and toroidal field has appeared in 
the region of strong toroidal magnetic field gradients near the top of the evolving toroidal 
magnetic layer. The instability appears in a quasi-two-dimensional manner with a wavevector 
primarily in the y—z plane , taking a form that is similar to a classical 'interchange' instability 



(ICattaneo h Hugheslll988l ). Roll-like motions principally swap lines of toroidal magnetic field 
that pierce y — z planes without a large degree of bending along the toroidal x-direction. 
There is a high degree of correlation between vertical velocity and toroidal field perturbations, 
strongly implying a buoyancy-driven instability (also see Figure [3]). 



We now compare C2 results with those of CI, which differs only the value of (. We 
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estimate the value of ( for C2 to be sufficiently large to render i?DD everywhere negative, 
and therefore anticipate that the instability found in CI, if double-diffusive, should not 
appear in this case. We run C2 for over twice the time that it takes for the instability in 
CI to manifest (up to t ~ 180) and the system simply continues to evolve according to the 
mean Equations (jSJ) and ([6]). Eventually, these dynamics are benignly disrupted by Alfvenic 
processes without any buoyancy instabilities (as explained in VB2) and we therefore halt 
the computation. We conclude that the instability discovered in CI is of a double-diffusive 
nature, since it depends critically on a sufficiently large thermal diffusivity (small £). 

A subtle issue in this problem is the influence of the strength of the imposed background 
field. As argued in VB1-2, this essentially provides a timescale for disruptive Alfvenic pro- 
cesses. For instability to occur, the necessary conditions must be met before Alfvenic pro- 
cesses can disrupt the source. In VB1-2, this required a strong velocity shear (hydrodynam- 
ically unstable, Ri ~ 0.03) in order to induce the required magnetic gradients sufficiently 
quickly for a regular (mean-density-deficit driven) magnetic buoyancy instability to occur. 
Things are more complex in the high Ri, double- diffusive regime considered here. To eluci- 
date this issue, we examine case C3 where the imposed background field strength is twice 
that of CI and C2. Figure ([3]) shows a vertical profile of the mean vertical transport of B x for 
both CI and C3, a good indicator of the existence and efficiency of any magnetic buoyancy 
instabilities. The measure plotted is 



n( \ wB x -wB x 

C ( Z ) = i R P ( 7 ) 

max z WrmsAi 



7 rms- L - , x,rms 



where the overline represents a horizontal average, and w TmB (z) = w 2 — w , B Xjims (z) = 
B 2 — B x . For ease of comparisons, the normalization factor in the denominator for all cases 
is that value calculated for CI and B x is rescaled to the same units. Figure ([3]) shows that 
the toroidal magnetic field perturbations correlate well with vertical velocity perturbations, 
confirming that the instabilities are most likely of a magnetic buoyancy type. However, the 
dotted curve in this figure shows that the transport is significantly weaker for C3 in the 
initial stages of the instability. This is a signature of the dynamic nature of the instability. 
Owing to the evolving background state, -Rdd does not entirely determine stability. In 
the equilibriu m stability problem for quasi-incompressible motion without an imposed shear 



considered in iHughes fc Weiss! (119951 ) . the growth rate of the instability becomes positive 
when -Rdd > 277r 4 /4. In the current problem, the portion of the domain that is unstable is 
perpetually moving as the magnetic layer evolves. Therefore, a locally unstable perturbation 
may not have enough time to amplify to a dynamically significant magnitude before it finds 
itself outside of the moving region of instability. Only if the growth rate of the instability 
is sufficiently large can an unstable mode amplify fast enough that the background appears 
effectively steady. Therefore, to exhibit instability in this evolving system unambiguously, 
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the growth rate of the instability must be at least of the order of the Alfvenic evolution 
timescale for the background, 

Pcrit (R T ,R B ,a,()~&z^/^(Q, (8) 

where, p rr \t (RT, Rr, cr, Q is the ma ximal growth rate as a function of the Rayleigh numbers 
as given in iHughes fc Weiss! (119951 ). and the Chandrasekhar number Q ^> 1. Equation (jSj) 
implies that in the presence of a stronger background field, the instability requires a faster 
growth rate to proceed in the same manner. Even though C3 has the exact same background 
B x {z) configuration as that for CI, the stronger background field strength, B Zj0 , alters the 
nature of the instability. At the location where CI grows strongly, C3 only reaches small 
amplitude (albeit in a shorter time: t « 44 cf. t m 88). However, the instability in C3 can 
still grow to significant amplitudes higher up in the domain as the system evolves further, as 
shown by the dashed curve in Figure ([3]). The amplitude is further modified since the stable 
background density stratification, and therefore i?DD changes slightly with height. 

We have demonstrated that the presence of shear-generated double-diffusive magnetic 
buoyancy instabilities can be controlled primarily by the single parameter (. Given the 
dynamic nature of the instability here, the second most important parameter is the Chan- 
dresekhar number, Q, which dictates the background Alfvenic timescales. C3 provides good 
evidence that the efficiency of the instability is governed through a relation at least qualita- 
tively similar to Equation ([8]). However, the question remains as to whether any other param- 
eters, such as the magnetic Prandtl number, ctm = cr/(, are critical to this process. Numerical 
constraints restrict our simulations to ctm — 0.5, which is less than unity thus preserving the 
correct orderi ng of timescales, but larger than the estimated solar value of % ~ 6 x 10~ 2 



(iGoughl 120071 ). To make some progress, we extend the analysis of IHughes h Weiss! ( 119951 ). 
noting again that this discusses the double-diffusive instabilites of an equilibrium rather than 
a dynamically-evolving system. It can be shown that Equation (|SJ) embodies the solution to 
a 12 th -order polynomial in p cvl t. In the limit 1 and R? -C this can be simplified to a 

cubic polynomial depending solely on the reduced variables r = Rt>t>/\Rt\, Q = Az 2 Q/\Rt\, 
and <jm'- 

AXl + X q X 2 r - 18X q X r - X q (27 + AX q ) = 0, (9) 

where 

x q = q ^ + ^\ x r = {r - q){1 + ^\ (10) 



Figure (j3J) shows the resulting dependence of the critical value of r versus q for several 
values of ctm- In particular, it can be seen that the instability is more easily obtained for 
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smaller <tm, which makes sense in terms of the lessening relative importance of viscosity, and, 
encouragingly, supports the viability of this process at more realistic solar parameters. 

The transport measure wB x — wB x has a more complicated dependence on o"m- In the 
appropriate Rayleigh number regime, wB x — WB X ~ p CTit + k% Tit , where k cr \ t is the horizon- 
tal wavenumber of the fastest growing mode. As <tm is decreased, k CT i t decreases and p CT i t 
increases, and hence the dependence on o"m is not straightforward. However, in the limit as 
cm ~~ ► 0, then wB x — wB x ~ p crit , which, from Figure PJ, is nonzero even at o"m = 0. This 
argument implies that, if we were able to run a computation at a significantly lower value 
of <7m than in CI and C3 (an extreme numerical challenge given our other parameters), we 
might expect that the buoyant transport should decrease, but not disappear altogether. 



Conclusions 



After the pessimistic results of VB1-2, it is encouraging that when double-diffusive 
effects are taken into account, the dynamics intuitively expected in the tachocline appear 
to be possible, satisfying the constraints that we have anticipated. We know roughly that 
( w 10~ 5 and Ri ~ 10 3 — 10 5 in the tachocline, and it is therefore plausible that their 
product is smaller than Az/H p m 1 (IGoughl 120071 ). However, we know very little about 
the configuration of any magnetic field in the solar tachocline or the associated Alfvenic 
timescales, and hence we do not know exactly how much smaller ( must be for the growth 
rate of an instability to become large enough for it to proceed. 

Even if instability occurs, it still remains unclear whether the nonlinear evolution of the 
system can produce magnetic fields that are sufficiently strong and coherent to be able to rise 
through the turbulent convection zone. While it now seems possible that high Ri shear flows 
could produce buoyant magnetic structures in the tachocline, one might expect that such 
buoyant structures would have magnetic energies that are (at most) in equipartition with the 
shear. If the turbulent convection zone acts locally in maintaining the tachocline's shear, then 
it is difficult to see how magnetic structures could obtain the significantl y greater strengths 



that are required to survive the disrupting effects of the convection zone (ITobias et al.ll2001 



Clinell2003l ; iFan et al.l 120031 ) . Of course, this question depends critically on how the energy 



of the shear relates to that of the convection zone, which is not well understood. 
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Table 1. Parameters for Simulations 



Case 


c 


a 


C K 


a 


Q 


Stability 


CI 


5.0 x 1CT 4 


2.5 x 10~ 4 


1.0 x 10~ 2 


1.25 x 10~ 5 


1.0 x 10 6 


unstable 


C2 


1.0 x 10~ 2 


5.0 x 10~ 3 


5.0 x 10~ 4 


1.25 x 10~ 5 


1.0 x 10 6 


stable 


C3 


5.0 x 10~ 4 


2.5 x 10~ 4 


1.0 x 10~ 2 


5.0 x 10~ 5 


4.0 x 10 6 


delayed onset 



Note. — In all cases, m = 1.6, 9 = 5, and Ri min = 2.96 
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Horizontally Averaged Toroidol Magnetic Field: B x (t,z) 
35 r ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — ■ — i — ■ — ■ — - 




Depth (Z) 



Fig. 1. — Time evolution of growing mean toroidal magnetic field vs. depth for CI and C2. 
The field is initially zero and builds a growing peak that ends at t « 88. 
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Fig. 2. — Volume-rendered images of flow and magnetic field from CI at t ~ 88. Image (a) 
shows vertical velocity, w (x, y, z) together with lines of magnetic field colored according to 
zonal velocity u(x, y, z) (red tones near the top of the box indicate flow of one direction, blue 
tones near the bottom of the box indicate opposite flow, and green tones indicate approximate 
stagnation). Image (b) shows fluctuating toroidal magnetic field, B x (x,y,z) — B x (z). The 
silver lines indicate streamlines of the fluctuating velocity (u — Uq, v, w). 
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Fig. 3. — Normalized buoyant vertical flux vs. depth for CI (solid) and C3 (dotted or 
dashed). All plots are proportional to wB x — wB x for each simulation and are scaled 
to be in the same units 



(B x rescaled, and all normalized by the same constant value 

(max 2 (w 2 — w 2 ^ (b 2 — B^ from CI). Since the Alfven timescales are twice as fast, 
the C3 plots are shown both when the instability is occurring in the same vertical location 
as CI (dotted, t p=* 44) and at the same actual time as the CI plot (dashed, t 88). 
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Fig. 4.— Critical R DD /\R T \ vs. Q/\R T \ for a M = a/( = 0,0.5,1.0,1.5,2.0, in the limits 
a, C < 1, and R T < 0. 



